%% The script for preprocessing the EEG data

% Author: Jia Liu
% Email: Jia.Liu@uts.edu.au

%%
eeglabDir = 'E:\Tool\eeglab14_1_2b\';
dataDir = 'C:\Users\EEG Data\';  % data folder
dataCleanDir = 'C:\Users\Processed\';  % data processed folder              

%% Create the subfolder to store the preprocessing datasets
if ~exist(dataCleanDir)
    mkdir(dataCleanDir);end

if ~exist(amicaDir)
    mkdir(amicaDir);end

%%
% General structure
% Step 1: Load the raw EEG dataset for each recording session. 
% Step 2: merge all sessions into one EEG. 
% Step 3: Downsample the data.            
% Step 4: High-pass filter the data at 1-Hz, and Low-pass fileter at 45-Hz. 
% Step 5: Apply clean_rawdata() to reject bad channels and correct continuous data using Artifact Subspace Reconstruction (ASR)
% Step 6: Run the ICA  
% Step 7: Estimate single equivalent current dipoles            
% Step 8: Save the cleaned dataset with ICA solution
% Step 9: Visualize the ICA solution

% Step 1: Load the raw EEG dataset for each recording session.
EEG1 = loadcurry([dataDir 'S01_1_active.cdt']);
EEG2 = loadcurry([dataDir 'S01_2_active.cdt']);
EEG3 = loadcurry([dataDir 'S01_3_active.cdt']);

% Step 2: merge all sessions into one EEG. 
EEG = pop_mergeset( EEG1, EEG2);
EEG = pop_mergeset( EEG, EEG3);

% Step 3: Downsample the data.
EEG = pop_resample( EEG, 250);
EEG = eeg_checkset( EEG );

% Step 4: High-pass filter the data at 1-Hz, and Low-pass fileter at 45-Hz.
EEG = pop_eegfiltnew(EEG, 'locutoff',1,'hicutoff',45,'plotfreqz',1);
EEG = eeg_checkset( EEG );

% Step 5: Apply clean_rawdata() to reject bad channels and correct continuous data using Artifact Subspace Reconstruction (ASR)
EEG = pop_clean_rawdata(EEG, 'FlatlineCriterion',5,'ChannelCriterion',0.8,'LineNoiseCriterion',4,'Highpass','off','BurstCriterion',20,'WindowCriterion',0.25,'BurstRejection','on','Distance','Euclidian','WindowCriterionTolerances',[-Inf 7] );
EEG = eeg_checkset( EEG );

% Step 6: Run the ICA  
runamica15(EEG.data, 'num_chans', EEG.nbchan,...
        'outdir', [dataCleanDir nameDataCleaned],...
        'num_models', 1, 'do_reject', 1, 'numrej', 15, 'rejsig', 3, 'rejint', 1);
EEG.etc.amica  = loadmodout15([dataCleanDir 'S01_active']);
EEG.icaweights = EEG.etc.amica.W;
EEG.icasphere  = EEG.etc.amica.S;
disp('runica successfull, storing weights and sphere.');
EEG = eeg_checkset(EEG, 'ica');
numICs = size(EEG.icaweights,1);

% Step 7: Estimate single equivalent current dipoles
[newlocs coordinateTransformParameters] = coregister(EEG.chanlocs, template_models(2).chanfile, 'warp', 'auto', 'manual', 'off');
templateChannelFilePath = [eeglabDir 'plugins/dipfit3.7/standard_BEM/elec/standard_1005.elc'];
hdmFilePath             = [eeglabDir 'plugins/dipfit3.7/standard_BEM/standard_vol.mat'];
EEG = pop_dipfit_settings( EEG, 'hdmfile', hdmFilePath, 'coordformat', 'MNI',...
	'mrifile', [eeglabDir 'plugins/dipfit3.7/standard_BEM/standard_mri.mat'],...
    'chanfile', templateChannelFilePath, 'coord_transform', coordinateTransformParameters,...
    'chansel', 1:EEG.nbchan);
EEG = pop_multifit(EEG, 1:EEG.nbchan,'threshold', 100, 'dipplot','off','plotopt',{'normlen' 'on'});
EEG = eeg_checkset(EEG);
    
% Step 8: Save the cleaned dataset with ICA solution
EEG = pop_saveset( EEG, 'filename', 'S01_active', 'filepath', dataCleanDir);
    
% Step 9: Visualize the ICA solution
pop_topoplot(EEG,0, [1:numICs] ,['S01_active_AMICA'],[] ,1,'electrodes','off');
print([dataCleanDir 'S01_active_AMICA.png'],'-dpng');
         